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ABSTRACT 

We study nonaxisymmetric perturbations of rotating relativistic stars, 
modeled as perfect-fluid equilibria. Instability to a mode with angular 
dependence exp(im0) sets in when the frequency of the mode vanishes. The 
locations of these zero-frequency modes along sequences of rotating stars 
are computed in the framework of general relativity. We consider models of 
uniformly rotating stars with polytropic equations of state, finding that the 
relativistic models are unstable to nonaxisymmetric modes at significantly 
smaller values of rotation than in the Newtonian limit. Most strikingly, the 
m=2 bar mode can become unstable even for soft polytropes of index N < 1.3, 
while in Newtonian theory it becomes unstable only for stiff polytropes of index 
iV < 0.808. If rapidly rotating neutron stars are formed by the accretion-induced 
collapse of white dwarfs, instability associated with these nonaxisymmetric, 
gravitational- wave driven modes may set an upper limit on neutron-star 
rotation. Consideration is restricted to perturbations that correspond to polar 
perturbations of a spherical star. A study of axial perturbations is in progress. 

Subject headings: ... 
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1. Introduction 

Rotating perfect-fluid equilibria are unstable to nonaxisymmetric instabilities driven 
by gravitational radiation. Although apparently damped by viscosity in neutron stars 
cold enough to have a superfluid interior, the instability may play a role in limiting the 
maximum rotation of newly formed neutron stars. In particular, if rapidly rotating neutron 
stars with weak magnetic fields form in the accretion-induced collapse of some white dwarfs 
(or the core-collapse of some stars), the instability may set an upper limit on rotation more 
stringent than the Kepler frequency - the frequency of a satellite in circular orbit at the 
star's equator. 

Prior to the present work, the onset of the nonaxisymmetric instability had been 
computed only in the Newtonian limit and estimated in the first post-Newtonian 
approximation and in a slow-rotation approximation.^ The first fully relativistic 
computation is presented here (and in Stergioulas' PhD thesis (1996)). 

That stars can be unstable to gravitational radiation was first found by Chandrasekhar 
(1970), who considered the m = 2 modes for Maclaurin spheroids (uniform density rotating 
stars) in a Newtonian context. Friedman & Schutz (1975, 1978a, 1978b), show that 
this instability also appears in compressible stars and that all rotating "stars" (rotating 
self-gravitating perfect-fluid configurations) are generically unstable to the emission of 
gravitational radiation. Even for slowly rotating models there will always be a polar mode 
of high enough mode number m (equivalently, of short enough wavelength) that is unstable. 

x After this was written, we received a preprint by Yoshida and Eriguchi (1997) 
computing instability points using a Cowling approximation in general relativity. This is an 
approximation that ignores changes in the gravitational potential, and it should be accurate 
for large-m modes. 
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For a perfect fluid, a nonaxisymmetric mode becomes unstable when its frequency vanishes 
in the inertial frame, i.e., with respect to an observer at infinity. This occurs when a star 
rotates fast enough that an inertial observer sees that a count errot at ing mode becomes 
corotating with the star. The amplitude of the perturbation then grows with time, making 
the star unstable. 

Realistic neutron stars are viscous, and the presence of viscosity will shift the onset 
of instability. Work by Lindblom & Mendell (1995) (see also Ipser & Lindblom 1991a,b; 
Yoshida & Eriguchi 1995) shows that viscosity will damp out the instability once the 
temperature is low enough that the interior is superfluid. A window apparently remains, of 
temperatures less than about 10 9 K and greater than about 10 10 K, where the instability 
will affect a sufficiently rapidly rotating neutron star as it cools down after its birth. In 
the Newtonian limit, it appears that only the I = m < 5 modes are not damped out by 
viscosity and the / = 4 mode apparently sets the most stringent limit on the maximum 
angular velocity; the critical angular velocity is about 90 % to 95 % of the Kepler limit. 
Old neutron stars, spun-up by accretion, may be too cold to be subject to the gravitational 
radiation driven instability. 

A surprise, recently pointed out by Andersson (Andersson 1997; Friedman and Morsink 
1997) is that axial modes for all values of m will be unstable for perfect-fluid models with 
arbitrarily slow rotation. In a spherical star, axial perturbations are time-independent 
convective currents that do not change the density and pressure of the star and do not 
couple to gravitational waves. For rotating stars, their growth time r is proportional to a 
high power of fl (rough arguments appear to imply r oc Q~ 4 ~ 2m ), and viscosity will again 
presumably enforce stability except for hot, rapidly rotating neutron stars. 

The onset of the nonaxisymmetric instability for several modes in rotating polytropes 
has been computed in the Newtonian limit by Imamura, Friedman & Durisen (1985) 
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using a Langrangian variational principle and by Managan (1985) and Ipser & Lindblom 
(1990) using an Eulerian variational principle constructed by Ipser & Managan (1985) . 
Cutler (1991) , Cutler & Lindblom (1992) and Lindblom (1995) estimate the correction 
to first post-Newtonian order. Weber et al (1991) provide an estimate in a relativistic, 
slow-rotation approximation. Since neutron stars are relativistic objects, relativity must 
have a significant effect on the onset of instability. The post-Newtonian analysis suggests 
that, for a given equation of state (EOS) and a given mode /, the ratio fl c /flk of the 
critical angular velocity where the mode becomes unstable to the maximum allowed angular 
velocity (the Kepler frequency Qk) decreases as the neutron star becomes more relativistic. 
Thus, nonaxisymmetric instabilities set a more stringent limit on the maximum angular 
velocity than Newtonian theory suggests. 

In this paper, we report the first computation of zero-frequency modes of rotating 
relativistic stars. Polytropes of index N = 1.0, 1.5 and 2.0 are considered, and attention 
is restricted to the fundamental (f) I — m < 5 modes, since these are the ones that are 
most likely to participate in limiting the maximum angular velocity of neutron stars. It is 
shown that the critical angular velocities of relativistic stars are considerably lower than the 
corresponding Newtonian estimates. As a result, the m = 2 bar mode can become unstable 
for much softer (N < 1.3) polytropes than in the Newtonian limit. This is expected to hold 
true for most realistic EOSs since they have an effective index N < 1.0. Thus, depending 
on the equation of state and the effect of viscosity, the m = 2 mode may participate 
in setting the upper limit on rotation for rapidly rotating neutron stars created by the 
accretion-induced collapse of white dwarfs. 

In forthcoming work we plan to include the effect of viscosity on the critical angular 
velocities and extend our method to iV < 1.0 polytropes and to realistic equations of state. 
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2. The Equilibrium Configurations 

The spacetime of a rotating relativistic star in equilibrium can be described by the 
stationary, axisymmetric metric of the form 



ds 2 = -e 2v dt 2 + e 2 ^{d<f) - udt) 2 + e 2a {dr 2 + r 2 d6 2 ), (1) 

(for a review see Friedman & Ipser (1992) ). 

In (HD, as in the rest of this paper, gravitational units (c = 1, G = 1) are employed. 
The metric involves four independent equilibrium metric potentials u, ip, a and u which 
are functions of r and 9 only. We assume uniform rotation, Q = constant, where Q is 
the angular velocity of the star. The matter is assumed to be a perfect fluid at zero 
temperature, described by a polytropic equation of state, for which the energy density e, 
pressure P and number density n satisfy the relations 



P = Kn 1+1 / N , (2) 
e = nm B + NKn 1+1/N , (3) 

where K is a constant, m-Q is the rest mass per baryon and iV is the polytropic index, 
related to the adiabatic index 7 by 

_d\nP_e + PdP_ 1 
7 ~ d\nn ~ P ~de ~ 1 + N' ^ 

The stress-energy tensor for a perfect fluid is 



T ab = (e + P)u a u b + Pg ab , 



(5) 
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where the equilibrium 4-velocity u a is given by 



u a = ^==(t a + nr), (6) 

V 1 — v 2 



with 



v = (n - u)e^~ v , (7) 

the fluid velocity with respect to a local zero- angular-momentum observer. We have 
denoted by t a and <p a the Killing vectors d t and d p hi associated with the time and rotational 
symmetries of the metric. 

An equilibrium model satisfies the field equation, 

Rab = 8vr (T ab - \^g a bT) , (8) 
and the implied equation of energy conservation 

u b V a T ab = 0^ u ^ a e = -( e + P)V a u a . (9) 

and Euler equation 

a bc X7 P 

q c b V a T ab = 0^ u "V a u b = ( ±-^, (10) 
Here q ab = g ab + u a u b is the projection of the metric tensor orthogonal to the 4-velocity. 

Numerical equilibrium models are constructed by a numerical code (Stergioulas & Friedman 
(1995) ) that implements the Cook, Shapiro & Teukolsky (1994) version of the KEH 
method (Komatsu, Eriguchi & Hachisu (1989) ). Our code has been shown to be accurate 
in an extensive comparison with other existing codes (Eriguchi et al. (1996) ). The field 
equations are solved on a 2-dimensional grid that is uniform in /i = cos 6* and in the radial 
coordinate 
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s = r/(r + r e ), (11) 

which compactifies the region r = to +oo to the finite region s = to 1 (r e is the 
value of the coordinate r at the equator). This is important, because the gravitational 
potentials have a nonnegligible value even far from a relativistic star; and it also simplifies 
the implementation of boundary conditions in the construction of both the equilibrium and 
the perturbed configurations. Note that the equator of the star is always at s = 0.5, so that 
(for a nonrotating star) half of the radial grid-points are inside the star and the other half 
are outside. 

We define dimensionless quantities by setting G — c — 1 and by fixing the length scale. 
If we define 



m B 

then K N I 2 can be used as the fundamental length scale. Dimensionless quantities will be 
denoted as e, Q etc. The dimensionless energy density and angular velocity are then 



K N G 
e = ; — e, 



(13) 



and 



K N / 2 

n = a (14 

c 

In subsequent sections, numerical results will be reported using dimensionless quantities 
as defined above. Note that the dimensionless quantities are independent of the polytropic 
constant K. 
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3. The Perturbed Configurations 

We will use an Eulerian formalism to describe fluid perturbations. The Eulerian 
perturbation in the metric tensor, 

5 g a b = Kb, , (15) 
can be determined by solving the perturbed field equations, 

5R ab = Sn(6T ab - ^g a bST - ^h ab T), (16) 

in a suitably chosen gauge. The perturbation in any other quantity is evaluated to first 
order in h ab . The change in the Ricci tensor is then (see e.g. Wald 1984 ) 

6R ab = V c V (Q /i c 6) - - (V a V b h c c + V c V c h ab ) , (17) 

where h a b = g ac h cb . Priou 1992 has explicitly computed the components of SR ab for a 
stationary, axisymmetric background in its most general form (prior to any choice of gauge). 
This allows us to directly use Priou's results with only the following modifications: 

i) since we are only interested in zero-frequency modes, we set the frequency of the mode 
and all time-derivatives equal to zero; 

ii) we rename Priou's functions, changing m to M, I to L and Q to y, in order to avoid 
confusion with the indices l,m that characterize a mode and with the angular velocity Q of 
the equilibrium star; and 

iii) we choose a gauge as described in section ||. 

The components of 5R ab in our gauge are displayed in Appendix |A]. 

The perturbation of the stress-energy tensor is 
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5T ab = u a u b (5e + 5P) + (e + P){u a 5u b + u b 5u a ) + g ab 5P + Ph ab , , (18) 
while the perturbation of its contraction is 

5T = 35P - 6e, (19) 

The change in pressure and energy density can be expressed in terms of a single scalar 
function SU, defined implicitly by the relation 

8P=(e + P)(8U + ^u a u b h ab ). (20) 
For an adiabatic perturbation, we have 

Se = { ^P, (21) 

and we will assume that the adiabatic index T of the perturbation be equal to the adiabatic 
index T of the equilibrium configuration, defined in (|j). 

In the perturbed relativistic Euler equations, 

S(q a c V b T bc ) = 0, (22) 

the only derivatives of Su a that occur are along the unperturbed fluid velocity u a . For 
perturbations with harmonic t- and ^-dependence, the equations are thus algebraic in 5u a 
and can be solved analytically (Ipser & Lindblom, 1992) . The perturbed fluid velocity is 
given by 



1 



8u a = % Q ab V b (5U + -u c u d h cd ) - 8F b + -u a u c u d h cd} (23) 



2 



1 



2 
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where 



SF a = q a b 



-u c V c (h bd u d ) + u c u d V b h cd + u c u d h cd 



e + P 



(24) 



The tensor Q ab involves only equilibrium quantities and the frequency a of the mode in a 
rotating frame; for the time-independent perturbations considered here, 



a = mVt. 



(25) 



Explicitly, 



Here 



Q ab = J_ (au t ) 2 q ab - 2uj a tt b + W(20 6 e ac fi c - <f> a e bc UJ, 



D 



t\2„ab o, ,o-r^b . ■ „,t In Ib^acc) la^bc, 



M * = M a V a t= \e 2v - e 2 ^(n - uf 



1/2 



(26) 



(27) 



the quantity 



D = (an*) 3 - 2au t Q a u a 



(28) 



is the determinant of Q 1 ; u a is the vorticity of the fluid, 



u* = e abcd u b V c u d ; 



(29) 



a is a generalization of the angular-velocity vector, 



n a = ^u t e abcd u b (v c t d + nv c <j )d y, 



(30) 



i a is the unit linear combination of the Killing fields that is orthogonal to u a , 



-q\<t> b ; 



(31) 
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and the tensor e is the volume element on the 2-surfaces orthogonal to the Killing 
trajectories 



e ab = e abcd <p c u d . (32) 

The tensor Q ab exists (and the perturbed Euler equations can be inverted) as long as 
the determinant in (|28T) does not vanish. 



For completeness, we list components of u a , u a and Q ab . With u l given by Eq. (|27|). 
we have 



u<p = Ua (j) a = e 2 ^(Q - w)u*; 



-if>—v—2a 



r(w*) 2 



„— tjj— v— la 



gtr = _grt = - iaf H>-v T J jy-^tf^ 

Q' e = -g e * = _j <7e -tf-« r - l u r D- x u% 

Q rr = D~ 1 [(au t ) 2 e~ 2a — (a/) 2 ] 

gr0 = _g0r = _ D -^J 

gr<A = _g^r = i af H>-'> r d>D-\l+Q.U t U+) 

Q ee = r -2 jD -l[( (JM *)2 e -2a_ r 2 ( ^ ) 2 ] 

ge0 = _q4,q = -i ae -*-» r -i u r D-i^i + Q u t u ^ 



The perturbed covariant fluid velocity is given by 
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5u a = Su b g ab + u h h ab . (33) 
Through (|20|) , (|2l|), (^3|), ( |24"1) and the equation (|3~3"D for the perturbed covariant fluid 



velocity, the perturbation in the stress-energy tensor ST ab is expressed entirely in terms of 
the perturbed metric h ab and the scalar 5U. Due to the lengthy substitutions involved, we 
used the algebraic program MAPLE to compute 5u a which was then substituted in 5T ab . 
The expressions for the components of ST ab in the stationary, axisymmetric background 
are several pages long (for each component), and we do not display them here. These 
expressions, along with the components of SR ab , are used to form the perturbed field 
equations (|16|). MAPLE is used to convert these expressions into numerical code for 
inclusion in an ANSI-C program. 

The perturbed field equations determine h ab for given 5U. The scalar function 5U is 
determined by an additional equation, the perturbed energy conservation equation 



5{u b V a T ab ) = 0. (34) 

which will be considered in section H. 



4. Expansion in Spherical Harmonics 

Because the equilibrium configuration is axisymmetric, linear perturbations of the star 
and geometry can be decomposed into a sum of terms with angular dependence e m< ^. We 
expect that, as is the case for spherical stars, each discrete mode of a nonrotating, spherical 
model has a continuous extension to a mode for each rotating model with the same equation 
of state. Modes of spherical stars have angular dependence given by the tensor, vector and 
scalar harmonics associated with a given YJ m (#, 0). The corresponding mode of a rotating 
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model can thus be labeled by I and m; but all harmonics Y™, with I' > I > m contribute to 
the mode. 

Studies of rotating Newtonian stars find nonaxisymmetric instability sets in first along 
a sequence with increasing rotation for a reflection-invariant polar mode with I = m, the 
lowest value of I for a given m. We have, in this first numerical study, correspondingly 
restricted our consideration to perturbations invariant under reflection in the equatorial 
plane (under the diffeo 9 — > n — 8)f\. 

The Eulerian perturbation in a reflection-symmetric, I = m, time-independent mode 
has the form 



rr oc Q>m+2l' Ym+2l' (35) 

l'=0 

oo 

h tr oc i a rn+2V+2 Y™ +2V+2 (36) 

l'=0 

oo 

Wfi OC Cl m +2l'+2 ^m+2l'+2 y, + &m+2/'+l $m+2Z'+l ft) (37) 

l'=0 

oo 

OC ^( a m+ 2V+2 ^m+2l'+2 ft + * &m+2/'+l $m+«'+l M ) (38) 
Z'=0 

oo 

^ OC ^ ( a m+2 «' $™ +2 Z' ,«/ + frm+2Z'+2 ^m+2l'+2 fiv 

+1 /ii/ 

(39) 



l'=0 



where the a's, 6's, and c's are real coefficients, different for each component of h a b- We have 
adopted the Regge- Wheeler (1957) notation for spherical harmonics: 



2 Recall that if a vector v a is invariant under the diffeo 9 — > 7r — its components 
f r = t> a V a r, f* = f a V a t, and = t> a V a are invariant, while t> e = t> a V a # changes sign 
because V a 9 changes sign. 
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e 


= d e Yr, 




— a r i ! 


(40) 


<&[ e 


= ~ sin 6^0^' 




= smOdeY™; 


(41) 


ee 






= 0; 


(42) 




= sin 2 6Y™; 






(43) 


Kee 


32vm 

— °e 1 \ ) 


iTrm 


= (dgd,/,- cot ^)17\ 


(44) 




= (<9|, + sin # cos ^YT"; 






(45) 


Xi ee 


= -L(^-cotW m , 


, ,m 

00 


4(sn!^ + C ° S ^ 


-sin fl<9 2 ) Y l m {AQ) 


Xi <$><$> 


= — sin ^(c^ — cot 6)d ( f ) Y l m . 






(47) 



In Priou's (1992) notation (and with our own redefinitions mentioned in section |3|), 
the perturbed metric h ab is expressed in terms of ten ^-dependent perturbation functions, 
h, p, k, w, q, a, b, L, M, and y, in the manner 



h tt 


= -2he 2u + (2yu + 2wcu 2 )e 2lp 


(48) 


h tr 


= L + aue 2 ^ 


(49) 


hte 


= M + bue 2 ^ 


(50) 


htcf. 


= -e 2lp (y + 2cuw) 


(51) 




= 2ke 2a 


(52) 


h r e 


= Q 


(53) 


h r( j) 


= -ae 2 ^ 


(54) 


hee 


= 2pr 2 e 2Q 


(55) 


h0<j) 


= -be 2 ^ 


(56) 




= 2we 2 ^, 


(57) 
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Equivalently, to first order in h a b 



ds 2 



e 2u (l + 2h)dt 2 + e 2 ^(l + 2w) d<f>-(u + y)dt - adr - bd6 



+e 



2a 



(l + 2k)dr 2 + r 2 (l + 2p)d9 2 



+ 2qdrdd + 2Ldtdr + 2Mdtd9. 



(58) 



An associated set of real, 0- independent, variables L, • • -w, can be defined by writing 



h 
k 
L 
M 

q 
y 

a 

b 

V 
w 



ke~ im4> , 
-i L e"™*, 
-i MsinO e- im<p , 
qsin6 e _<m * 
ysm 2 9 e- im<t >, 
-i asin 2 e _im * 
-i bsm 3 6 e _im * 
psm 2 6 e- im<t> , 
wsin 2 e e- im<t> . 



(59) 
(60) 
(61) 
(62) 
(63) 
(64) 
(65) 
(66) 
(67) 
(68) 



Then, the metric perturbations of these modes can be expanded as follows in terms of 
Legendre polynomials: 



h, k, L, y, a, p, w ~ sin m ^ a 2 i>(r)P 2 i> (cos 0), 

l'=0 



(69) 



-17- 



M,q,b ~ sin m # £a 2r+1 (r)P 2r+1 (cos£), (70) 



l'=0 

It is obvious from (|69|) and (70|), that the functions M, g and 6 will be antisymmetric 
about the equatorial plane, while all other functions will be symmetric. Furthermore, with 
these definitions, the time-independent perturbed field equations become a system of ten real 
equations for ten real unknowns. 

The corresponding expansions for axial I = m modes, for which the tensor h a b changes 
sign under reflection (modes that reduce to axial modes in the Regge- Wheeler gauge in the 
spherical limit), can easily be written in a similar fashion. We note that the behavior under 
reflection in the equatorial plane will be opposite to that of polar modes. 



5. Gauge Choice 

A linear perturbation of a star, described by the set of quantities (5e, dp, 5u a , h a b) 
is physically equivalent to the gauge-related perturbation described by the set (8e + C v e, 
Sp + C v p, 5u a + C v u a , h ab + £ v g a b) for any smooth vector field rf that preserves the 
asymptotic behavior of the metric. The choice of gauge is important for the successful 
numerical solution of the perturbed field equations. After experimenting with a large 
number of possible gauges, we found that a numerical solution was more easily obtained in 
a gauge defined by the four conditions 



h r g 
hoc/, 
ht(/> 



= 

= 

—u)h, 
he. 



q = 0, 
b = 0, 



y = o, 



2{ip-a) 



W = p. 



(71) 
(72) 
(73) 
(74) 



-18- 



The components of h a t, satisfy all imposed boundary conditions in this gauge. In 
addition, we required that it reduce to the Regge- Wheeler polar gauge in the non-rotating 
limit, in which only h tt , h rr , h$Q and are non-zero, with Kqq and satisfying the 
equivalent of ( f74|) in a Schwarzschild metric (see previous section). With this choice of 
gauge there are six nonzero metric functions in the list fl6l])-(|68|), namely h, p, k, L, M 
and a. They will be expressed in terms of the function 5U by solving six components of 
the field equation 5R a b = 8n5(T a b — ^g a bT), (tt), (rr), (99) (tr), (t9) and (r9). Note that 
condition (|74|) implies that the perturbation function p does not have to be redefined, as in 
(|67|), because it has the desired angular behavior 

oo 

p ~ sin m # 5> 2I ,(r)P 2P (cos0). (75) 

l'=0 

In Appendix ^| we list the necessary components of the perturbation in the Ricci tensor 
in this gauge. 



6. Numerical Solution 

The task of solving the coupled system of six differential equations is not trivial. In the 
gauge specified in section || the (tt) and (99) equations are elliptic for h and p respectively. 
The (rr) equation is parabolic for k (it is missing a d 2 k/dr 2 derivative). The (tr) and (r<j)) 
equations are second order ordinary differential equations (ODEs) for L and a in the angular 
direction, while the (t9) equation is a second order ODE for M in the radial direction. Each 
type of equation requires its own finite-difference scheme and boundary conditions. Still, 
we were able to solve all six equations simultaneously on a 2-dimensional finite grid. 

In the radial direction, the grid coincides with that of the equilibrium configuration; 
that is, it consists of a number of grid-points, uniformly spaced in the coordinate s, which 
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ranges from s = (r = 0) to s = 1 (r = oo), the equator always being at s = 0.5 
(the coordinate s is defined by (|TTD). In the angular direction a different grid is used. 
Following Ipser & Lindblom (1990) , the n angular grid-points are located at the angles 
(ii = cos Qi which correspond to the zeros of the Legendre polynomial of order 2n — 1: 
-P2n-i(/ i i) = 0. This has the advantage of using only a small number of angular grid-points 
to describe the star, which results in a smaller linear system of equations to be solved. Since 
the perturbation functions are either symmetric or antisymmetric when reflected in the 
equatorial plane, one need only solve for < 9 < n/2, and the boundary conditions at 
9 = 7r/2 are incorporated in the expressions for the angular derivatives. A drawback of this 
method is that stars with stiff equations of state (N < 1.0) are less accurately described 
than stars with soft EOSs (N > 1.0), because the former have discontinuous derivatives of 
energy density and metric functions across the surface and this results in the appearance of 
Gibbs phenomena in the angular derivatives. 

The above choice of angular grid-points allows us to use high-order formulae for the 
angular derivatives. In all six equations, all first and second order angular derivatives are 
approximated as 



-ll~f ± (s i ,» j ) = it D tkf ± (si,»k), (76) 
" fc=i 



and 



^fHsi^j) = J2 H fkf ± (s i ^k), (77) 

°f J ' A=l 



where f^ 1 is a function which is either symmetric (+) or antisymmetric (— ) under the 
transformation 9 — ► n — 9, i.e. / ± (s,— jj) = ±/(s, fi). The functions Dp k and Hp. are 
derived in Appendix p] and are constructed specifically for functions that have the angular 
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behavior 



/+ = sin m 9 y £a !a ,P 2l r(jA), 

l'=0 



(7i 



like h, k,p, L and a, and 



f- = sin m 0$> 2m P 2m (//), (79) 

l'=0 

like M. 

In all equations, except for the parabolic equation (rr) , the first and second order radial 
derivatives are approximated by the standard, second order accurate, central difference 
formulae 



If. 

ds 1 '-' 



fi+l,j fi- 

2As 



1,3 



10) 



and 



d r _ fi+lj - 2fi,j + fi-l,j / R1 \ 

where As is the distance between two radial grid-points. In the parabolic equation (rr), 
we use an implicit solution scheme; that is, we use a first order accurate, backwards, first 
radial derivative for k, 



dS f ^ = As ' (82) 

while for the first order derivatives of other functions involved in the (rr) equation we 
still use the second order accurate equation (|80D . Second order radial derivatives in the 



(rr) equation are not approximated by flST]) but by a similar equation using twice the 
grid-spacing 



d 2 r fi+2,j - 2fi,j + fi-2,j 

ds* h > 3 (2As) 2 ' 1 } 



If one tries to use (q1|), the solution for k oscillates. That using (|^) instead of (|8l|) can 
suppress such oscillations was discovered in the construction of equilibrium models, where 
a second order radial derivative appears in the source term of a non-elliptic equation, by 
Stergioulas & Friedman (1995) . Mixed derivatives involved in some of the equations are 
approximated by 



d 2 1 n 

d^ f ^ = A7£^ (Afc ~ (84) 

At the center and at infinity the appropriate boundary finite difference formulae are used 
in all equations. 

For the two elliptic equations (rr) and (88) we require that the solution has vanishing 
first order angular derivative at 8 = tt/2 and vanishes at infinity. The boundary condition 
at 8 = 7r/2 is built into the construction of the angular derivative formulae. No boundary 
conditions are needed at the center and on the symmetry axis for the elliptic equations, but 
the discretized equations force the solution to vanish there for any mode with m ^ 0, as is 
the case for the exact solution. 

For the two second order angular ODE's, (tr) and (r0), the boundary conditions at 
8 = 7r/2 and on the symmetry axis are set by the angular derivative equations ( |T0D and 
(\T1\). For the second order radial ODE (t8) we require that M vanishes at the center and 
at infinity. Finally, for the parabolic equation (rr), boundary conditions at 8 = n/2 and on 
the symmetry axis are set by the angular derivative equations ([F6|) and ([f?]), while in the 
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radial direction we set k = at the center only (at infinity the boundary is open). 

Of the above equations, the most difficult to solve numerically are the (rr) and (r<f>) 
equations. In the (rr) equation, the solution is propagated from the center to infinity via 
a first order radial derivative, so it is far more sensitive to local inaccuracies (like those 
produced by the finite grid-spacing at the surface) than the elliptic equations. In particular, 
the surface of the rotating star does not correspond to a constant value of the radial 
coordinate but falls in between grid points. For polytropes of index N > 2, the energy 
density goes smoothly enough to zero at the surface that this does not pose any problem 
even with a grid of only 200 radial points. For polytropes of index 1 < iV < 2, however, a 
small jump in the solution for k appears at the surface. Stiff polytropes of index N < 1 
have, in addition, discontinuous second order derivatives of the equilibrium metric functions 
across the surface and a jump in the numerical solution for k is unavoidable when using our 
grid, since the angular derivative equations were designed only for smooth functions. The 
small jump appearing in the numerical solution for k does not affect any other perturbation 
function significantly, except for the function a. Another problem with the (rr) equation - 
of same origin - is that, for stiff equations of state, the solution oscillates in the vacuum 
region, unless one makes the approximation 



The error introduced by making this approximation is small, since for all stars that we 
examined 



d 2 h 




d 2 h 



(85) 



ds 2 



ds 2 



h ~ —p. 



(86) 



Thus, the significant benefit of making this approximation (suppression of all oscillations in 
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the solution) outweighs its mild cost. 

Finally, the (r<f>) equation seems to be very sensitive to local inaccuracies, especially at 
the surface of the star and near its center. The inaccuracies near the center of the star occur, 
because the perturbation function a depends on the differences between other perturbation 
functions, which all have very small values (compared to their maximum value) near the 
center. However, the perturbation function a is of lesser importance compared to the other 
five functions, when computing the critical angular velocities. In fact, we determined that a 
only minimally affects the other five perturbation functions, so that setting it equal to zero, 
results in almost unchanged solutions for the other perturbation functions and unchanged 
critical angular velocities for even the most relativistic stars of the stiffest polytropic index 
(N = 1.0) we examined. 

We also determined that the other two off-diagonal perturbation functions L and M 
are weakly coupled to the diagonal ones. This is expected, since in our gauge the perturbed 
metric reduces to a diagonal form in the non-rotating limit and is the dominant 
off-diagonal perturbation (g t( p is the only non-vanishing off-diagonal metric element for a 
rotating star). We emphasize that all the above statements hold for perturbations of the 
type considered in the present paper. 

The finite-differencing of the system of perturbed field equations yields a large linear 
system. The unknowns are ordered in a way that casts the matrix of the linear system 
in a band-diagonal form. A direct solution method (LU-decomposition) is employed. 
Alternatively, we also use an iterative Bi-Conjugate Gradient method with a symmetric 
succesive over-relaxation (SSOR) method as a preconditioner. The direct solution was 
faster, but required a larger amount of random access memory.]] 

3 The subroutines for the solution of the linear system we used, are part of the Portable 
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A representative solution of the coupled system of six perturbation equations, for a 
specific trial function SU, is shown in Figs. (HD - ©■ The background star is a rapidly 
rotating, N=1.5 polytrope with dimensionless central energy density e c = 0.061 and 
dimensionless angular velocity Q = 0.092. The solution was obtained for the I = m = 3 
mode using a trial function 

5U = r l Pr(v) } (87) 

where P™(/x) is the associated Legendre function. This trial function is the dominant 
term in the true eigenfunction of the I = m modes (see section |8]). All six equations were 
solved on the same 101 x 7 (radial x angular) grid. The perturbation functions h, p, k, L, 
M and a are shown along all seven angular spokes and at the pole (where they vanish). 
With the exception of M, the perturbation functions are symmetric under reflection in the 
equatorial plane and their maximum value occurs for 6 = tt/2. The perturbation function 
M is antisymmetric under reflection, so it vanishes at 9 = tt/2. Its maximum value occurs 
about half-way between the pole and the equator. The diagonal perturbation functions h, p 
and k have no nodes inside the star. This is a characteristic of ha for fundamental I = m 
modes in the Newtonian limit. The perturbation function k (the solution to the parabolic 
(rr) equation) exhibits a small jump at the surface, as was explained earlier. In addition, 
near the surface, the solution can have a wave-like character. This is normal and arises 
because the solution is obtained for a trial function SU and not for the true eigenfunction. 
The off-diagonal functions L, M and a vanish much faster in the vacuum surrounding the 
star than the diagonal functions do; they become negligible almost immediately outside the 
surface of the star. The function L does not have any node inside the star, but M and a 

Extensible Toolkit for Scientific Computation (PETSc) package developed at Argonne 
National Laboratory (Gropp and Smith 1994) . 
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do. Finally, the solution for the function a exhibits some oscillations at the surface of the 
star, which are induced by the behavior of the function k there. Regarding other polytropic 
indices and modes, the perturbation functions remain similar to those in Figs. @) - (§)• 
For stiffer polytropes and for more relativistic stars, the solutions are peaked closer to the 
surface than for softer polytropes and less relativistic stars. For larger m, the solutions are 
more narrowly peaked about their maximum value because the dominant radial behavior of 
SU is r l . For N = 1.0 polytropes the oscillations at the surface are larger than for N = 1.5 
polytropes, while for N = 2.0 polytropes they are negligible. Otherwise, the solutions for 
the perturbation functions are very similar to those in Figs. (P - (0). 

7. A Truncated Gauge 

In the Newtonian limit, the important component of the perturbation is h a , which 
reduces to — 2h. In addition, the metric functions h, k, p and w satisfy the relation 



For energy densities typical in neutron stars, the h a component still dominates h^. We 
determined numerically, that even for the most relativistic and rapidly rotating N — 1.0 
polytrope, the relation (|88|) still holds approximately in our gauge. This led us to investigate 
whether by making certain approximations, one can still obtain accurate values for the 
critical angular velocities of neutral modes, while solving fewer than six equations. Indeed, 
we find that by making the approximations 



h = —k = —p = —w. 




h tt h, 



rr 



h = —k = —p, 



(89) 




g tt + 2ug t< f, g. 



rr 



9ee 



and 
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h te = h r4> = M = a = 0, (90) 

one can obtain very accurate critical angular velocities, while solving only two perturbed 
field equations. Equation ( |g9D is motivated by (|5Ej) and fl9"0|) is used since t^g = = 
for the equilibrium configuration and because h t e, h r $ vanish in the spherical limit in the 
Regge- Wheeler polar gauge , i.e. they are proportional to the background metric function 
u). In this truncated gauge, the perturbation in the metric tensor, can be written in the form 



( -2h(e 2u - u 2 e 2 *) L -2oohe 2 ^ \ 

-2he 2a 

sym. —2hr 2 e 2a 

V -2he 2 ^ J 



(91) 



Only two perturbation functions are nonzero, h and L, which are determined by solving the 
{tt) and (tr) equations. 

Table [l] shows a comparison of critical configurations obtained in the full and truncated 
gauges. The equilibrium star is a N = 1.0 polytrope with e c = 0.3, which is close to the 
central energy density of the maximum mass configuration. The critical configurations for 
the I — m — 3, 4 and 5 modes are computed in both gauges on the same 101 x 7 grid. 
We compare the dimensionless critical angular velocity and the critical ratio of rotational 
to gravitational binding energy T/|VK| C in the two gauges. The critical values of T/|VK| 
differ by less than 1% and the critical angular velocities by less than 0.5% when computed 
in the truncated gauge compared to computing them in the full gauge. The advantage 
of using the truncated gauge is that far less memory is needed by the numerical code. 
While six equations can be solved in a reasonable time with a maximum 201 x 12 grid 
on a DEC Alpha with 256 MBytes of memory, two equations can be solved with a much 
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finer grid of 801 x 12 or more points. The execution speed in not limited by the processor 
speed (the matrix inversion is very fast) but by the available random access memory For 
very soft polytropes of index iV > 2.0, a 201 x 12 grid gives sufficiently accurate critical 
angular velocities, so that one can use the full gauge for these models. But, for polytropes 
ofiV<1.0a201xl2 grid determines the critical angular velocity with an error that can 
exceed several percent or more. This was determined by comparing critical configurations 
for Newtonian polytropes, obtained with our code, to published results. Thus, for realistic 
neutron stars, it is necessary to obtain the critical angular velocities with a grid of at least 
801 x 12 points and this could be achieved only in the truncated gauge. For this reason, all 
results reported in the present paper were obtained in the truncated gauge with a grid of 
801 x 12 points, unless otherwise stated. 



8. The Perturbed Energy Conservation Equation 

In preceding sections it was shown how a solution of the perturbed field equations 
is obtained for a given trial function SU. The complete description of the neutral modes 
requires the satisfaction of both the perturbed field equations ([16]) and the perturbed energy 
conservation equation ([34]) by the perturbation in the metric h a b and the scalar function 
SU. This system of equations fll6|) and ( j34|) has a zero-frequency solution only for the stars 
for which a zero-frequency mode exists. For any equilibrium model, however, one can solve 
only the perturbed field equations (16) for trial functions SU and use the perturbed energy 



conservation equation fl34 ) to construct a criterion for locating the marginally stable star 



(for a given mode) along a sequence of rotating stars. 

The perturbed field equations fllB"]) are implicitly linear in the function SU, as is the 
perturbed energy conservation equation (|34|). Hence ( ]34|) can be represented as a linear 
operator L acting on a function SU 
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L(5U) = 0. (92) 



Equation (92) represents an eigenvalue problem for the linear operator L with zero 



eigenvalue and eigenfunction SU. 

In the Newtonian limit, it was shown that the eigenfunction SU of an I = m neutral 
mode can be approximated accurately by expanding it in terms of a set of basis functions 
SUi, 

SU = J2a i 5U i , (93) 

i 

of the form 

SUi = SU^ k) = r l+2 ^Y™ 2k (cos 9), (94) 

with j, k = 0, 1, ... and with each set of indices (J, k) yielding a particular SU^. In practice 
the eigenfunction SU is represented with reasonable accuracy by only a few terms, and the 
(0,0) term r z YJ m dominates the expansion (pl|) . 



Substituting the expansion (|93|) in ( ^2]) yields 

a-MSUi) = 0. (95) 

i 

If we define the inner product^ 



4 The motivation for defining the inner product as in (|96D will become apparent in 
Appendix p. 
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< 6Uj\L\6Ui > = J i-^L{8Ui)y/^d 3 x, (96) 

where g is the determinant of the equilibrium metric tensor, then taking the inner product 
of ( p5|) with respect to SUj yields 

^2a i <6U j \L\6U i >=0. (97) 

i 

Although the basis functions are not orthogonal to each other, non-zero < 5Uj\5Ui > terms 
do not appear on the r.h.s. of ( p7| ) because we require that the eigensystem has zero 
eigenvalue. The last equation is a linear, homogeneous system for the coefficients c^. The 
system has a non-trivial solution only if its determinant vanishes. This yields a criterion for 
locating the zero-frequency modes 

Criterion A stationary, axisymmetric model has a nonaxisymmetric, zero-frequency mode 
with angular dependence e im ^ , if 

det < 5Uj\L\8Ui >= 0. (98) 



In practice, we start with a slowly rotating star and compute the matrix elements 
< 5Uj\L\5Ui >. For slowly rotating stars the determinant in (|98"D always has a large value. 
Keeping the central energy density constant, we look at stars of increasing angular velocity, 
until the determinant goes through zero. The star for which (|98|) is satisfied is the one for 
which the particular I = m mode has a zero-frequency solution and the nonaxisymmetric 
instability sets in through that mode. In this method, the accuracy in locating the neutral 
modes depends on how well the expansion fl9"3| ) approximates the true eigenfunction SU. 
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In Appendix y, we obtain, in terms of Eulerian quantities, an expression for 

< SUj\L\SUi > that is used in the numerical computations. As discussed in this Appendix 
the above method of solving the perturbed energy conservation equation is not identical 
to using a variational principle; because of our choice of field equations, the matrix 

< 5Uj\L\5Ui > is not symmetric. For central energy densities typical in neutron stars, 
however, the matrix is nearly symmetric and the method nearly coincides with a variational 
principle. 

9. Critical configurations 

Following the method developed in previous sections, we computed the neutral mode 
(critical) configurations for the fundamental I = m = 2, 3, 4 and 5 nonaxisymmetric modes. 
Three equations of state were examined, having polytropic indices N = 1.0, 1.5 and 2.0 
polytropes. The N = 1.0 polytropes include models with mass and radius similar to those of 
realistic neutron stars. The other two equations of state were examined for completeness and 
for comparison with published results in the Newtonian limit. The critical configurations 
were computed in the truncated gauge with a fine grid of up to 801 x 12 (radial x angular) 
grid points. The trial function 6U was expanded as in (p4[) using different number of terms. 
Eight basis functions, corresponding to the indices j = 0, ...,3 and k = 0, 1 in (|94|) were 
sufficient to determine the critical configurations with good accuracy, especially for the 
iV = 1.5 and iV = 2.0 polytropes, where the error consistently decreases with increasing 
number of basis functions and increasing number of grid points. This was not so for the 
N — 1.0 polytropes. Owing to the finite grid-size, the surface of the star is discontinuous, 
and this gives rise to Gibbs phenomena at the surface when using the expansions (|76|), 
([T?|) for the angular derivatives. As a result, increasing the number of grid-points did not 
monotonically decrease the error. This behavior is expected for our choice of coordinates 
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and grid. It is the price one has to pay for being able to solve a smaller than otherwise 
linear system of equations. A similar behavior is reported in Bonazzola et al. (1996) , who 
use a similar angular grid in computing the neutral, viscous bar mode in relativistic stars. 
Still, by computing the N — 1.0 configurations with various grid sizes and various numbers 
of basis functions, we could obtain sufficiently accurate results. 

Our code was checked in the Newtonian limit by comparing the critical configurations 
for the / = m = 3, 4 and 5 modes for N = 1.0, 1.5 and 2.0 polytropes to results published 
by Managan (1985) , Imamura, Friedman & Durisen (1985) and Ipser & Lindblom (1990) 
. As can be seen in Table [2] we are in very good agreement with all published Newtonian 
results, the exception being T/|W| C for the m = 5 mode in N = 1.0 polytropes, which 
differs by 2 % from the value published in Ipser & Lindblom (1990) , from which Imamura, 
Friedman & Durisen (1985) differ by 3 % (in the opposite direction). This highlights the 
increased inaccuracy involved in computing neutral modes for N < 1.0 polytropes. 

For relativistic polytropes we obtain the following results: For N = 1.0, figure |] shows 
the ratio of the critical angular velocity Q c to the Keplerian angular velocity Qk at same 
central energy density as a function of central energy density for the four modes examined. 
The lowest central energy density in the figure corresponds to a mildly relativistic star. The 
highest central energy density shown is the central energy density of the most massive (and 
thus most relativistic) star allowed by the particular equation of state. The filled circles on 
the left vertical axes represent the values of Q c /Qk in the Newtonian limit. As the central 
energy density increases and the star becomes more relativistic, Q c /Qk decreases and it 
decreases at a faster rate as it approaches the most relativistic configuration. Contrary to 
the Newtonian limit, where N = 1.0 polytropes do not have an unstable m = 2 mode, 
we find that in relativistic N = 1.0 polytropes the m = 2 mode becomes unstable when 
the central energy density exceeds roughly one 10th the central energy density of the most 
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massive star. At the most relativistic configuration, the m = 2 mode becomes unstable 
for T/|iy | c =0.065 or Q c /Qk — 0.91. Table ^ compares the most relativistic critical 
configurations to their counterparts in the Newtonian limit. The value of Q c /Qk decreases 
by roughly 15% for the m = 3, 4 and 5 modes. Figure (|5|) shows the critical ratio T/|W| C 
for the same N — 1.0 polytropes. This ratio decreases faster and by a larger percentage 
than Q c /Qx, owing to the fact that the Keplerian value T/\W\k also decreases as one 
samples more relativistic stars. In Table |3| it is shown that for the most relativistic N = 1.0 
polytrope the critical ratio T/|W| C is about 40 % smaller, for the m = 3, 4 and 5 modes, 
than the corresponding ratio in the Newtonian limit. When the decrease in the Keplerian 
value Ty|W|jr is taken into account and one looks at the ratio (T/\W\ C )/ (T/\W\k) the 
most relativistic values are still 25 — 30% lower than the Newtonian values. 

Figures (§) and (0) and Table ^ display the critical configurations for N = 1.5 
polytropes. These polytropes are softer than the ones with index N = 1.0. Consequently, 
they have a smaller maximum mass and are less relativistic. For this reason, relativity has 
a smaller effect on the onset of the nonaxisymmetric instability. The m = 2 mode does 
not become unstable even for the most relativistic N = 1.5 polytropes. For the m = 3, 4 
and 5 modes, the value of Q c /Qk decreases by 7 — 10% for the most relativistic models 
compared to the Newtonian limit. The corresponding decrease for T/|W| C is 30 — 35% and 
for (T/\W\ C )/(T/\W\ K ) it is 13 - 19%. 

Plots of Q c /Qk and T/\W\ C for the N = 2.0 polytropes are shown in figures (|8p and 
These are extremely soft models and their maximum mass occurs at nearly Newtonian 
central energy densities. Again, the m = 2 mode does not become unstable. As seen in 
table g, for the m = 3, 4 and 5 modes the values of Q C /Q K , T/\W\ C and (T/\W\ C )/(T/\W\ K ) 
are only a few percent less than in the Newtonian limit, for the largest central energy 
density. Polytropes constructed with iV = 2.0 do not resemble realistic neutron stars, but 
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are included here for completeness. 

Regarding the accuracy of our results, we estimate that the determined critical angular 
velocities and critical T/\W\ C are accurate to better than 2 % for N = 1.0, 1-2 % for 
N = 1.5 and 1 % for N = 2.0 polytropes. 

10. Discussion 

We have treated nonaxisymmetric neutral modes in the context of general relativity. 
We have found a gauge in which six coupled perturbed field equations can be solved 
simultaneously with good accuracy. Furthermore, we found an approximate gauge, in which 
the critical configurations for iV > 1.0 polytropes are located with sufficient accuracy, 
while solving only two independent perturbed field equations. We showed that general 
relativity has a large effect on the location of nonaxisymmetric neutral modes and forces 
the nonaxisymmetric instability to set in for smaller rotation rates than Newtonian theory 
suggests. 

The large effect of relativity on the onset of the nonaxisymmetric instability is most 
striking in the case of the m = 2 modes. In the Newtonian context it was shown that 
uniformly or nearly uniformly rotating neutron stars cannot become unstable to the 
m = 2 mode unless the equation of state is excessively stiff (see e.g. Skinner & Lindblom 
(1996) ). For polytropes, the classical result by James (1964) restricts the onset of the 
m=2 instability to polytropes having an adiabatic index larger than r crit = 2.237 (which 
corresponds to a polytropic index A^ crit = 0.808). We find that, in the context of general 
relativity, this critical value becomes r crit = 1.77 (polytropic index A^ cri t = 1.3). 

In the Newtonian limit, the m = 2 mode driven by gravitational radiation coincides 
with the m=2 mode driven by viscosity even for compressible fluids, as was shown by Ipser 
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& Managan (1985) . It is interesting that work by Bonazzola et al. (1996) suggests that 
in general relativity, the m = 2 viscosity-driven mode has a critical adiabatic index only 
slightly higher than James's result. Thus, in general relativity, the viscosity driven and the 
gravitational radiation driven m = 2 modes no longer coincide and the effect of relativity 
seems to be very different on each of them. 

In forthcoming work, we plan to study realistic equations of state and include the effect 
of viscosity on the onset of the nonaxisymmetric instability. 
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very helpful discussions. This research has been supported by NSF grants PHY-9105935 
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A. Perturbed Ricci Tensor 

We list the perturbation in the six components of the Ricci tensor that are used in 
the perturbed field equations, in the gauge of section [|. We follow the notation introduced 
in section [|. In the r.h.s. of each equation, subscripts denote partial derivatives, e.g. 
h H = d 2 h/d(j) 2 . 
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2 [ 2r z 1 

+ ^e^- 2 ^^(c 2 + ^) + |^e 2 ^)(4^ — ktf, — h+) 

1 2W-a-,)K (2a0 + ^ _ 3 ^ )M _ 2 L + ^ / _ 3 ^ _ 1 

2 [ 2 \ r 

+^(L 9 w e -M^ r )-Iw rr -^M| (A6) 



B. Angular Derivative Formulae 

In this appendix we derive high-order, finite-difference formulae that approximate the 
angular derivatives of functions that are known at the discrete set of angles fii = cos 9i 
(with i = l...n), which correspond to the zero's of the Legendre polynomial P2 n -i{^i) = 0. 
A function f(r,fi)e im< ^ can be expanded in terms of associated Legendre polynomials as 

00 

/(r )A i) = E/*(0^M(A*). ( B1 ) 

k=0 

Since PJ^X\ m \ is a polynomial of order k multiplied by (1 — /i 2 )'" 1 '/ 2 , expansion (pl|) is 
equivalent to 

00 

/(r >A i) = (1-// 2 ) H/2 E /!(0n(A*), (B2) 

k=0 



-38- 



where Pk{fj) is the Legendre polynomial of order k. By orthogonality of the Legendre 
polynomials, the coefficients /|(r) in the expansion ( |B2|) are determined as 



fl(r) = (k + -) jj{T^){l-^y^PMd^. (B3) 
The integral in ( (FJ3|) can be computed with high accuracy using Gaussian quadrature 

/ 9 + (r, fi)dfi ~ J2 w i9 + ( r , Hi), ( B4 ) 
where g + (r,fi) is a function symmetric in fi and are weights that are tabulated in e.g. 



Abramowitz & Stegun . For a function g (r, /i), antisymmetric in /z, the integral in (|B4j ) 
vanishes. 

We are interested in functions f{r,fi) that have definite reflection symmetry, i.e. that 
are of the form f ± (r,fi) = ±/ ± (r, — fi). For f + (r,fi), the nonvanishing coefficients in (|B3|) 
become 



n i 

4 = E~( 4fc + 1M1 - tfr lm]/2 P2k(Vi)f + (r, m), (B5) 



while for / (r, //), flB"3|) yields 



n i 

4 + i = E 2^ 4fc + 3 )^( X - ^ 2 ) _H/2 ^+i(^)/-(r,/x,). (B6) 



i=l 



Differentiating (B2) with respect to fi one obtains 



d_ 



\m\fM 



d 



+ (i-^) H/2 E/iM^W, 



(B7) 



fc=0 
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and substituting ( P5|) into ( |B7|) yields 



|"/ + (r,^) = EA5/ + (r,M,), (B8) 



where 



D± 



-5 i:j + (1 - ^ 



2\|m|/2 



71-1 I Q 

x: -(4*: + - ^ 2 )- |mi/2 p 2fc (/i J )^p 2 ,(/i i ), 



(B9) 



and 5ij is the Kronecker delta. Next, we use the recurrence relation 







2k 



in (|B9| ) and define 



(BIO) 



n—l 



k=l 



(Bll) 



In (|B9] ) we truncated the expansion to include Legendre polynomials up to order P 2 n-2 (/•*)• 
If we repeat this for f~(r, fx) and define 



n ~ 2 (2k + 1) 

S ij = J2 o ^ + 3 ) jP 2fc+l(A i i) -HiP2k+l(Hi) + P2k(^i) 



fc=0 



then D^j is given by 



13 (1 " A# 



\m\ni5ij + 



\m\/2 



WjSfj , . 



and the first-order angular derivative formula becomes 

l:/ ± (r,A*i) = f:i^/ ± (r,^). 

j=i 



(B12) 



(B13) 



(B14) 
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In order to obtain the second-order angular derivatives, we differentiate (B8) with 
respect to \i 



dfx- 



\m\(l + nf) + m 2 /i 2 

oo 

+ (i-ft) w/ TA 



+ 2\m\tM d + 



l-//f <9/T 



d 2 



Differentiating the recurrence relation ( |Bf U| ) with respect to fi we obtain 



(B15) 



d 2 2k r r i 

P2k{^) = - — 2 { -[l + (l-2A0^]P 2 *0O 



dfx 



(1 - ^: 



+ (3 - 4A;) y u J P 2fe _ 1 (/i i ) + (2k - l)P 2fe _ 2 (^) }. 



Continuing in the same fashion as for the first-order derivatives, we define 



(B16) 



ij (i - i4Y 



1 + (1 + |m|)// 2 \m\8. 



'ij 



-2|m|^(l-/z 2 )P± + 



1*1*1/2 



(B17) 



where 



n-l 



Gg = £A;(4A;+l)P 2fc (^) { -[l + (l-2k)fi 2 ]P 2k (fi t ) 

+ (3 - 4A;)/i i P 2fe _ 1 ( y u i ) + (2* - l)P 2k - 2 (Hi) }, 



k=i 



(B18) 
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and 



n-2 



^ = E (2fc 9 +1) ( 4fc + 3 ) P 2 fc+ i(^) { -(l-2ktf)P 2k+1 ( f i i 



k=0 



2 



+ (1 - 4/c)^P 2fe (^) + 2A;P 2fc _i }. 

(B19) 

Then, the second-order angular derivatives for the functions / ± (r, /i) are 

1 f±(r,^) = Y, H tf ± (r,H)- (B20) 



<9/i 2 



C. Perturbed Energy Conservation Equation 

In this Appendix, we construct a variational principle for the perturbed energy 
conservation equation (based on Friedman & Ipser (1992) ) and show its relation to the 
method presented in section [8]. In addition, we derive an explicit expression for the inner 
product < 5Uj\L\5Ui > defined in section ||. 



C.l. A Variational Principle 

In the Lagrangian formalism, perturbations are described by the Lagrangian 
displacement vector £ a . The complete description of a perturbation requires the solution of 
the perturbed field, Euler and energy conservation equations for the metric perturbation 
h a b and the displacement vector £ a . If one can solve all but one of the above equations, 
using for example a trial vector £ a , then the remaining equation will not be satisfied. One 
then has to construct a criterion, which will be used to identify the equilibrium star for 
which the remaining equation is satisfied. 

For perturbations (£ a , h a b) that satisfy the perturbed energy conservation equation 



5(u b V a T ab ) = 0, 



(CI) 



the perturbed field equations form a symmetric system (Friedman & Schutz (1975) ). That 
is, two pairs (£ a , h ab ) and (£ a , h ab ) satisfying ( |CT| ) obey the symmetry relation 

i b 5(v c T bc ) + ^-h bc 5{G bc - 8nT bc ) = -C(i a , h ab ; C, h ab )+V b R b , (C2) 

where G ab = R ab — ^g ab R c c and V b R b is a divergence constructed from £ a , h ab and 
h ab . In ( |C2|) , C(^ a , h ab ; h ab j is a function symmetric under the interchange of the trial 
solutions (£ a , h ab ) and (<f a , h ab )\ hence the r.h.s. in (|C2] ) is symmetric up to a divergence. 

In the Eulerian approach to solving the perturbation equations, the perturbed Euler 
equation 



S(q a b V a T bc ) = 0, (C3) 

is solved analytically and six components of the perturbed field equation are solved for h ab 
given a trial function 5U. The symmetry in ( p2[) can be exploited in order to construct 
a variational principle for the remaining unsolved equations. The solved perturbed Euler 
equation can be eliminated from (|C2] ) by decomposing the Lagrangian displacement vector 
£ a into vectors normal and parallel to the 4-velocity 



£ a = G. - &u c )u a , (C4) 
where = Q a b£, b - Equation (|C3|) then implies 
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i b 6(v a T ab ) = -{i c u c ) 5(u b V a T ab ). (C5) 
For trial solutions the perturbed energy conservation equation is not satisfied 

5(u b V a T ab ) = L(5U) ± 0, (C6) 

where L is a linear operator acting on the function SU. 

The Lagrangian displacement £ a has a gauge freedom in its component along u a : 
Adding a vector field fu a (where / is some arbitrary scalar function) to £ a leaves the 
Eulerian perturbations unchanged. As was shown by Friedman & Ipser (1992) , if (|C1| ) is 
not satisfied, the perturbation equations can still be cast in a symmetric form as in ( |C2| ) 
(but with a redefined divergence term V b R b ) if the component of £ a along u a is given the 
value 



1 



-5U. (C7) 
iau l 

With this definition, 



4 6(v a T ab ) = l -^L{5U). (Q 

(711 



Next, we define 



F{5U; SU) 



1 

167T 
1 

16/T 



''be 



8ttT 



be 



be 



9 bc h d d ) 5 



Rbr — 



87i(T bc - -g bc T d a 



(C9) 



-44- 



which is implicitly bilinear in SU and SU. The gauge freedom in h ab leaves only six 
(out of ten) components of the perturbed field equations independent. Thus, if a trial 
solution satisfies six components of the perturbed field equations in (|C9|) , the other four 
components of the perturbed field equations will be an implicit functional of the perturbed 
energy conservation equation, the only remaining equation that needs to be satisfied. 
Schematically, one can write 



T(SU\ SU) = SUF L(5U) 



(CIO) 



where F is a functional of L(SU). We have used the linearity of ^(Stl; SU) to factor out 
SU. The symmetry relation ( C2|) becomes 



— C(£ a , h ab ; £ a , h ab 



iSU 



SU 



L(SU) + SUF L(SU) 



V h R b 



— - + F L{SU) - V b R b 



= -C(SU;SU) 



(Cll) 



A variational principle for the perturbed energy conservation equation is constructed, 
by requiring that the following integral (which is implicitly quadratic in SU) vanishes 



1 = J -C(SU; SU)^gd 3 x = 0. 
Because the integral of the divergence vanishes, this yields 



(C12) 







i = Jsu 







L(SU)^/=gd 3 x = 0. 

The integral in ( |C13| ) is stationary with respect to first order variations in SU 



(C13) 
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51 

5{5U) 



aw 
L(6U) = 0, 



L(SU) = 0, 



(C14) 

(C15) 
(C16) 



provided that ^L(5U) ^ — F L(SU) . Thus, / = is a variational principle for the 
perturbed energy conservation equation L(5U) = 0. 



In practice, one can expand the trial function 5U in terms of a set of basis functions 



SUi 



5U = J2 aM. 



Substituting flCl7D into flC13l) yields 



(C17) 



]T J2 a i a 3 I 5U j \— t + F L(5Ui)^^d?x = 



<7U 



a i a i A {a) = °> 



* 3 



where we defined the symmetric matrix A with elements 



(C18) 



A i3 = I 8U 3 



+ F 



aw- 



L(5Ui)^gd z x. 



(C19) 



Equation ( |C14| ) implies that the integral I will be stationary to the variation of any of the 
coefficients Oj in the expansion ( U17[ ). Thus, 



5I_ 

5di 



J2 a j A (v) = °- 



(C20) 
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The last equation is a homogeneous linear system for the coefficients aj which has a 
nontrivial solution only when 

detA {ij) = 0. (C21) 

Since we can find an explicit form for Ay in terms of known quantities, we could have 
used (|(J21|) as a criterion for locating neutral modes. However, the term F[L(8Ui)] involves 



many second order angular and radial derivatives and it is not certain how accurate its 
evaluation on our finite grid would be. Since the method described in section S gives a 
substantially simpler criterion which does not involve F[L(5Ui)}, we chose to use that for 
locating the neutral modes. It is interesting to note that the matrix elements 



<6U i \L\6U i >= ( l -^-L{5U l )^jd z x (C22) 

J <7U 

used in section |8| are nearly symmetric under the interchange of SUi and 8Uj, for all 
configurations considered. This indicates that the method described in section |8| nearly 
coincides with a variational principle. 

C.2. An Expression for the Inner Product 

Since for the equilibrium star u b V a T ab = 0, 



L(SU) 



5[u b V a T a 
A{u b V a T ab 



-A 



[e + P)V b u b + u b V b e 



-u c V c [Ae + -(e + P)q ab Ag ab 
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Ae+~(e + P)q ab Ag ab 



(C23) 



(cf. Friedman & Ipser, 1992). Then, 



L(6U) = -iau 1 { 5e + £* V a e + ^ (e + P) h c c + u a u b h ab + 2 V a f a + 2uV V & } , (C24) 



where we used Ag ab = h ab + V a £& + Vb£ a . Using the decomposition (|C4[) of £ a , one obtains 



(C25) 



and 



«VV a 6 = SU - C±u a V a u b . 



(C26) 



The Euler equations for the equilibrium configuration yield 



u V b u c 



qJ^bP 
' e + P ' 



(C27) 



Then (|C23| ) becomes 



L(6U) = -iauH 5e + eiV a (e + P) + ~(e + P) h c c + u a u b h ah + 2V 



(C28) 



The perturbed energy density is 



= ( e + P ) 2 ( tfc/ + 



(C29) 
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so that 



L(5U) = -iauH (e + P) 



{e + P). TT If, e + P\ a h c c 

-6U + -1 H hiV/i a6 + — 

pr 2 1 PT 2 



+ V a £i(e+P) • (C30) 



The matrix elements defined in section Bl now take the form 



< SUjlLlSUi > 



J{e + P)Uu i 



(j+g) 
PT 



PT 



(C31) 



where we have used the time-independence of £ a to eliminate the term 
/ ^a[£i( e + P)5U]y/^ ~gd 3 x as an integral of a spatial divergence. 

Finally, the component of £ a normal to the 4-velocity u a is related to the component of 
5u a normal to u a by 



8. 



5u a , 5u a - \u b u c h hc u a 



IOU 1 



iau l 



(C32) 



(cf. Ipser & Lindblom (1992) ) and the expression for < 5Uj\L\5Ui > used in our numerical 
computations becomes 



<8U j \L\8U i > = J{t + P)\5Uj 



^ 5U ^\{ l + t -pf 



(5u a - \u b u c h bc u a ) 



law 



P\ h c 



-gd 3 x, 



(C33) 



where h ab and 5u a are computed with SUi. 
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Table 1: Comparison of critical configurations in the full and truncated gauges 

N = 1.0 

e c = 0.3 T/\W\ C Q c 

7 x 101 



Full gauge 4.63e-2 2.82e-l 

Truncated gauge 4.59e-2 2.81e-l 



m=4 



Full gauge 3.46e-2 2.48e-l 

Truncated gauge 3.47e-2 2.48e-l 



m=5 



Full gauge 2.83e-2 2.26e-l 

Truncated gauge 2.83e-2 2.26e-l 
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Table 2: Comparison of critical T/|W| in Newtonian limit with other authors 







m=3 


m=4 


m=5 




N= 


1.0 






Present 




7.92e-2 


5.79e-2 


4.62e-2 


Managan (1985) 




7.94e-2 


5.81e-2 




Imamura et al. (1985) 




8.0e-2 


5.8e-2 


4.4e-2 


Ipser & Lindblom (1990) 


8.00e-2 


5.84e-2 


4.53e-2 




N= 


1.5 






Present 




5.61e-2 


4.33e-2 


3.36e-2 


Managan (1985) 




5.6e-2 


4.3e-2 




Imamura et al. (1985) 




5.7e-2 


4.3e-2 






N= 


:2.0 






Present 




3.35e-2 


2.81e-2 


2.28e-2 


Managan (1985) 




3.3e-2 


2.8e-2 
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Table 3: Critical configurations for N=1.0 polytropes 



e c t/\w\ c (t/\w\ c )/ n c n c /n K 

(T/\W\ K ) 



m=2 



Relativistic 3.4e-l 6.49e-2 0.777 3.44e-l 0.911 



m=3 



Newtonian 1.0e-8 7.92e-2 0.769 6.69e-5 0.921 
Relativistic 3.4e-l 4.55e-2 0.544 2.96e-l 0.783 



m=4 



Newtonian 1.0e-8 5.79e-2 0.562 5.94e-5 0.818 
Relativistic 3.4e-l 3.53e-2 0.422 2.64e-l 0.699 



m=5 



Newtonian 1.0e-8 4.62e-2 0.449 5.41e-5 0.745 
Relativistic 3.4e-l 2.87e-2 0.343 2.40e-l 0.635 
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Table 4: Critical configurations for N=1.5 polytropes 



e c T/|W| C {T/\W\ C )/ Q c Q C /Q K 
(T/\W\ K ) 



m=3 



Newtonian 1.0e-7 5.61e-2 0.943 1.62e-4 0.980 
Relativistic 6.1e-2 3.82e-2 0.804 1.02e-l 0.917 



m=4 



Newtonian 1.0e-7 4.33e-2 0.728 1.47e-4 0.886 
Relativistic 6.1e-2 2.79e-2 0.587 8.87e-2 0.800 



m=5 



Newtonian 1.0e-7 3.36e-2 0.565 1.32e-4 0.796 
Relativistic 6.1e-2 2.34e-2 0.492 8.18e-2 0.738 
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Table 5: Critical configurations for N=2.0 polytropes 



e c T/\W\ C {T/\W\ C )/ Q c Q C /Q K 
(T/\W\ K ) 



m=3 



Newtonian 1.0e-8 3.35e-2 0.991 3.67e-5 0.997 
Relativistic 5.1e-3 2.59e-2 0.951 2.22e-2 0.979 



m=4 



Newtonian 1.0e-8 2.81e-2 0.832 3.42e-5 0.928 
Relativistic 5.1e-3 2.10e-2 0.771 2.02e-2 0.895 



m=5 



Newtonian 1.0e-8 2.28e-2 0.676 3.12e-5 0.848 
Relativistic 5.1e-3 1.70e-2 0.624 1.84e-2 0.814 
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Fig. 1. — Representative solution for the perturbation functions h and p (at = 0). Notice 
that p ~ —h even for this very relativistic configuration. The scaling of the vertical axis is 
determined by the trial function 8U. The equator of the star is at s = 0.5, while s = 1.0 
corresponds to infinity. The solution in each figure is shown at pi = cos# = 0, 0.23, 0.45, 
0.64, 0.80, 0.92, 0.98 and 1.0. The maximum is at pi = 0. 

Fig. 2. — Representative solution for the perturbation functions k (at (f) — 0) and a . Notice 
that k ~ —h. The scaling of the vertical axis is determined by the trial function SU. The 
equator of the star is at s = 0.5, while s = 1.0 corresponds to infinity. The solution in each 
figure is shown at pi = cos 9 = 0, 0.23, 0.45, 0.64, 0.80, 0.92, 0.98 and 1.0. The maximum is 
at pi — 0. 

Fig. 3. — Representative solution for the perturbation functions L and M. The scaling of 
the vertical axis is determined by the trial function SU. The equator of the star is at s = 0.5, 
while s = 1.0 corresponds to infinity. The solution in each figure is shown at pi = cos# = 0, 
0.23, 0.45, 0.64, 0.80, 0.92, 0.98 and 1.0. For L the maximum is at pi = 0. For M the dashed 
line is at pi — 0.23 and the maximum is at pi — 0.45. 

Fig. 4. — Critical angular velocity over Keplerian angular velocity at same central energy 
density vs. the dimensionless central energy density e c for the m = 2, 3, 4 and 5 neutral modes 
of N = 1.0 polytropes. The largest value of e c shown corresponds to the most relativistic 
stable configurations, while the lowest e c corresponds to less relativistic configurations. The 
filled circles on the vertical axis represent the Newtonian limit. 

Fig. 5. — Critical ratio of rotational to gravitational binding energy vs. the dimensionless 
central energy density e c for the m = 2, 3, 4 and 5 neutral modes of N — 1.0 polytropes. The 
largest value of e c shown corresponds to the most relativistic stable configurations, while the 
lowest e c corresponds to less relativistic configurations. The filled circles on the vertical axis 
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represent the Newtonian limit while the dotted line is the Kepler limit. 

Fig. 6. — Critical angular velocity over Keplerian angular velocity at same central energy 
density vs. the dimensionless central energy density e c for the m = 3, 4 and 5 neutral modes 
of N — 1.5 polytropes. The largest value of e c shown corresponds to the most relativistic 
stable configurations, while the lowest e c corresponds to less relativistic configurations. The 
filled circles on the vertical axis represent the Newtonian limit. 

Fig. 7. — Critical ratio of rotational to gravitational binding energy vs. the dimensionless 
central energy density e c for the m = 3, 4 and 5 neutral modes of N = 1.5 polytropes. The 
largest value of e c shown corresponds to the most relativistic stable configurations, while the 
lowest e c corresponds to less relativistic configurations. The filled circles on the vertical axis 
represent the Newtonian limit while the dotted line is the Kepler limit. 

Fig. 8. — Critical angular velocity over Keplerian angular velocity at same central energy 
density vs. the dimensionless central energy density e c for the m = 3, 4 and 5 neutral modes 
of N = 2.0 polytropes. The largest value of e c shown corresponds to the most relativistic 
stable configurations, while the lowest e c corresponds to less relativistic configurations. The 
filled circles on the vertical axis represent the Newtonian limit. 

Fig. 9. — Critical ratio of rotational to gravitational binding energy vs. the dimensionless 
central energy density e c for the m = 3, 4 and 5 neutral modes of N = 2.0 polytropes. The 
largest value of e c shown corresponds to the most relativistic stable configurations, while the 
lowest e c corresponds to less relativistic configurations. The filled circles on the vertical axis 
represent the Newtonian limit while the dotted line is the Kepler limit. 
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